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Abstract 

Previous  analyses  of  composite  clocks  have  been  presented  with  2-state  clock  models  for 
random  clock  phase  and  frequency  deviations.  But  Kalman  filter  estimate  errors  due  to  3-state 
clock  frequency-drift  deviations  have  not  been  previously  presented.  I  have  employed  the  3- 
state  Zucca-Tavella  clock  model  to  simulate  an  ensemble  of  four  3-state  clocks.  I  have  found  a 
dominant  common  component  for  3-state  clock  Kalman  filter  phase  errors  with  the  curvature  of 
a  nonlinear  low-order  polynomial  that  does  not  exist  for  2-state  clocks.  Cesium  clocks  are  free 
of  frequency-drift,  but  rubidium  and  hydrogen-maser  clocks  have  significant  frequency  drift. 
Use  of  the  3-state  clock  model,  and  an  understanding  of  its  true  and  estimated  behavior,  will 
facilitate  operation  of  the  associated  composite  clock. 


I.  INTRODUCTION 

GPS  Time  is  created  by  processing  GPS  pseudo-range  measurements  with  the  operational  GPS  Kalman  filter. 
Brown  [2]  refers  to  the  object  created  by  the  Kalman  filter  as  the  GPS  Composite  Clock,  and  to  GPS  Time  as  the 
Implicit  Ensemble  Mean  phase  of  the  GPS  Composite  Clock.  The  fundamental  goal  by  the  USAF  and  the  USNO 
is  to  control  GPS  Time  to  within  a  specified  bound  of  UTC/TAI.  I  present  here  a  quantitative  analysis  of  a 
simulated  GPS  Composite  Clock,  derived  from  detailed  simulations  and  associated  graphics.  GPS  clock  diffusion 
coefficient  values  used  here  were  selected  to  enable  a  characterization  of  3-state  clocks  that  illuminate  the  effects 
of  frequency-drift.  Diffusion  coefficient  values  used  for  each  of  the  four  clocks  S 1 ,  S2,  N 1 ,  and  N2  are  presented 
in  several  of  the  initial  figures.  SI  and  S2  refer  to  ground  station  clocks.  N1  and  N2  refer  to  NAVSTAR  clocks. 

My  interest  in  the  GPS  Composite  Clock  derives  from  my  interest  in  performing  real-time  orbit  determination1  for 
GPS  NAVSTAR  spacecraft  from  ground  receiver  pseudo-range  measurements.  The  estimation  of  NAVSTAR 
orbits  would  be  incomplete  without  the  simultaneous  estimation  of  GPS  clock  parameters.  I  use  simulated  GPS 
clock  phase,  frequency,  and  frequency-drift  deviations,  and  simulated  GPS  pseudo-range  measurements,  to  study 
Kalman  filter  estimation  errors. 

I  am  indebted  to  Charles  Greenhall  (JPL)  for  encouragement  and  help  in  this  work. 


II.  COMPLETE  ESTIMATION  AND  CONTROL  PROBLEM 

The  USNO  operates  two  UTC/TAI  master  clocks,  each  of  which  provides  access  to  an  estimate  of  UTC/TAI  in 
real-time  (1  pps).  One  of  these  clocks  is  maintained  at  the  USNO,  and  the  other  is  maintained  at  Shriever  Air 


371 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

NOV  2007 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2007  to  00-00-2007 

4.  TITLE  AND  SUBTITLE 

Composite  Clocks  with  3-State  Models 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Analytical  Graphics,  Inc, 220  Valley  Creek  Blvd, Exton, PA, 19341 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

39th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting,  26-29  Nov  2007,  Long  Beach,  CA 


14.  ABSTRACT 

Previous  analyses  of  composite  clocks  have  been  presented  with  2-state  clock  models  for  random  clock 
phase  and  frequency  deviations.  But  Kalman  filter  estimate  errors  due  to  3-state  clock  frequency-drift 
deviations  have  not  been  previously  presented.  I  have  employed  the  3-state  Zucca-Tavella  clock  model  to 
simulate  an  ensemble  of  four  3-state  clocks.  I  have  found  a  dominant  common  component  for  3-state  clock 
Kalman  filter  phase  errors  with  the  curvature  of  a  nonlinear  low-order  polynomial  that  does  not  exist  for 
2-state  clocks.  Cesium  clocks  are  free  of  frequency-drift,  but  rubidium  and  hydrogen-maser  clocks  have 
significant  frequency  drift.  Use  of  the  3-state  clock  model,  and  an  understanding  of  its  true  and  estimated 
behavior,  will  facilitate  operation  of  the  associated  composite  clock. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

20 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


39th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


Force  Base  in  Colorado  Springs.  This  enables  the  USNO  to  compare  UTC/TAI  to  the  phase  of  each  GPS  orbital 
NAVSTAR  clock  via  GPS  pseudo-range  measurements,  by  embedding  a  UTC/TAI  master  clock  in  a  USNO  GPS 
ground  receiver.  Each  GPS  clock  is  a  member  of  (internal  to)  the  GPS  ensemble  of  clocks,  but  the  USNO  master 
clock  is  external  to  the  GPS  ensemble  of  clocks.  Because  of  this,  the  difference  between  UTC/TAI  and  the  phase 
of  each  NAVSTAR  GPS  clock  is  observable.  This  difference  can  be  (and  is)  estimated  and  quantified.  The  RMS 
(Root  Mean  Square)  on  these  differences  quantifies  the  difference  between  UTC/TAI  and  GPS  Time.  Inspection 
of  the  differences  between  UTC/TAI  and  the  phase  of  each  NAVSTAR  GPS  clock  enables  the  USNO  to  identify 
GPS  clocks  that  require  particular  frequency-rate  control  corrections.  Use  of  this  knowledge  enables  the  USAF  to 
adjust  frequency  rates  of  selected  GPS  clocks.  Currently,  the  USAF  uses  an  automated  bang-bang  controller3  on 
frequency-rate. 


III.  STOCHASTIC  CLOCK  PHYSICS 

The  most  significant  stochastic  clock  physics  are  understood  in  terms  of  Wiener  processes  and  their  integrals. 
Clock  physics  are  characterized  by  particular  values  of  clock-dependent  diffusion  coefficients,  and  are 
conveniently  studied  with  aid  of  a  relevant  clock  model  that  relates  diffusion  coefficient  values  to  their  underlying 
Wiener  processes.  For  my  presentation  here,  I  have  selected  “The  Clock  Model  and  Its  Relationship  with  the 
Allan  and  Related  Variances”  presented  as  an  IEEE  paper  by  Zucca  and  Tavella  [18]  in  2005.  Except  for  FM 
flicker  noise,  this  model  captures  the  most  significant  physics  for  all  GPS  clocks.  I  simulate  and  validate  GPS 
pseudo-range  measurements  using  simulated  phase  deviations,  simulated  frequency  deviations,  and  simulated 
frequency-drift  deviations  according  to  Zucca  and  Tavella.  Figure  1  presents  an  ensemble  of  simulated  random  3- 
state  clock  phase  deviations,  overlaid  with  independently  calculated  3ct  boundaries. 


Clock  Phase  Simulation  Validation 


Zucca-Tavella 


Figure  1.  3-State  clock  phase  simulations. 
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IV.  GRAPHICS  FOR  CASE  13 

True  deviations  for  3-state  clocks  were  simulated  for  Case  13  and  were  used  to  construct  GPS  pseudo-range 
measurements.  The  measurements  were  processed  by  the  Kalman  filter  (KF1)  described  below  with  initial  state 
estimate  clock  deviation  errors.  The  simulated  (true)  clock  phase  deviations  and  estimated  clock  phase  deviations 
are  presented  in  Figure  2.  Usually,  one  expects  the  phase  estimates  to  lie  close  to  the  truth.  This  is  not  the  case 
here,  because  the  clock  phase  deviations  are  not  observable  from  the  GPS  pseudo-range  measurements. 
Nonetheless,  estimated  clock  phase  deviations  are  created  by  KF1.  The  simulated  (true)  clock  frequency 
deviations  and  estimated  clock  frequency  deviations  are  presented  in  Figure  3.  Inspection  of  Figure  3  helps  to 
explain  the  divergence  with  time  between  true  clock  phase  deviations  and  estimated  clock  phase  deviations 
presented  in  Figure  2. 


Simulated  &  Estimated  Phase  Deviates  (s) 


Case  13 


Figure  2.  Simulated  &  estimated  phase  deviates. 


The  estimated  clock  deviations  were  subtracted  from  the  true  (simulated)  clock  deviations  to  enable  inspection  of 
KF1  estimation  errors.  Estimation  errors  for  the  four  clocks  in  phase,  frequency,  and  frequency  drift  are  presented 
in  Figures  4,  5,  6,  and  7.  Figure  4  is  deceptive  in  that  one  might  be  led  to  believe  the  KF1  phase  errors  are  smooth 
with  time.  Figure  5  magnifies  an  interval  of  Figure  4  to  demonstrate  that  KF1  phase  errors  are  not  smooth  for  3- 
state  clocks.  Except  for  initial  KF1  phase  errors,  Figure  4  demonstrates  that  most  of  the  phase  estimation  errors 
for  the  four  clocks  SI,  S2,  Nl,  and  N2  are  similar.  These  similar  phase  errors  are  the  unobservable  components 
of  KF1  estimation  error  common  to  each  clock.  Figure  5  enables  identification  of  KF1  estimation  errors  that  are 
independent  for  each  clock  with  perturbations.  These  perturbations  are  the  observable  components  of  KF1 
estimation  errors.  Continued  processing  of  GPS  pseudo-range  measurements  will  reduce  the  variances  on  the 
observable  components. 
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The  common  unobservable  phase  component  for  the  four  clocks  SI,  S2,  Nl,  and  N2  has  been  estimated  with  a 
second  Kalman  filter  (KF2).  This  common  phase  component  was  subtracted  from  the  KF1  phase  error  for  SI  to 
enable  identification  of  the  observable  component  in  KF1  phase  error  for  SI.  The  latter  is  presented  in  Figure  8. 


V.  KALMAN  FILTER 

I  present  my  approach  for  the  optimal  sequential  estimation  of  clock  deviation  states  and  their  error  covariance 
functions.  Sequential  state  estimates  are  generated  recursively  from  two  multidimensional  stochastic  update 
functions,  the  time  update  (TU),  and  the  measurement  update  (MU).  The  TU  moves  the  state  estimate  and 
covariance  forward  with  time,  accumulating  integrals  of  random  clock  deviation  process  noise  in  the  covariance. 
The  MU  is  performed  at  a  fixed  measurement  time  where  the  state  estimate  and  covariance  are  corrected  with  new 
observation  information. 

The  sequential  estimation  of  GPS  clock  deviations  requires  the  development  of  a  linear  TU  and  nonlinear  MU. 
The  nonlinear  MU  must  be  linearized  locally  to  enable  application  of  the  linear  Kalman  MU.  Kalman’s  MU  [9] 
derives  from  Sherman's  Theorem  [12,13,10];  Sherman's  Theorem  derives  from  Anderson's  Theorem  [i];  and 
Anderson's  Theorem  derives  from  the  Brunn-Minkowki  Inequality  Theorem  [14,5].  The  theoretical  foundation 
for  my  linearized  MU  derives  from  these  theorems. 


Simulated  &  Estimated  Frequency  Deviates  (s/s) 
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Figure  3.  Simulated  &  estimated  frequency  deviates. 
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VA.  Initial  Conditions 

Initialization  of  all  sequential  estimators  requires  the  use  of  an  initial  state  estimate  column  matrix  X(l  0  and  an 
initial  state  estimate  error  covariance  matrix  P(^Q  for  time  tQ . 

VB.  Linear  TU  and  Linear  MU 

Derivation  and  calculation  for  the  discrete-time  Kalman  filter,  linear  in  both  TU  and  MU,  is  best  presented  by 
Meditch  [10],  Chapter  5. 


KF1  Phase  Errors  &  KF2  UECC  Estimate  (s) 


Case  13 


Figure  4.  KF1  phase  errors  &  KF2  UECC  estimate. 


VC.  Linear  TU  and  Nonlinear  MU 

The  simultaneous  sequential  estimation  of  GPS  clock  phase  and  frequency  deviation  parameters  can  be  studied 
with  the  development  of  a  linear  TU  and  nonlinear  MU  for  the  clock  state  estimate  subset.  This  is  useful  to  study 
clock  parameter  estimation,  as  demonstrated  in  Section  VIII. 
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Let  X  ■i  denote  an  n  *  1  column  matrix  of  state  estimate  components,  where  the  left  subscript  j  denotes  state 
epoch  tj  and  the  right  subscript  i  denotes  time-tag  tt  for  the  last  observation  processed,  where  i  ,  j  e  {0,1,2,...}. 
Let  Pjt  denote  an  associated  n  x  n  square  symmetric  state  estimate  error  covariance  matrix  (positive 
eigenvalues). 


KF1  Phase  Errors  &  KF2  UECC  Estimate  (s) 


Figure  5.  KF1  phase  errors  &  KF2  UECC  estimate  magnified. 


Linear  TU 

For  k  e  {0,1,2,. .  ,,M},  the  propagation  of  the  true  unknown  n  x  1  matrix  state  X  k  is  given  by: 

X k+\  =  ^k+\,k^k  +  Jk+\,k 

where  Jk+l  k  is  called  the  process  noise  matrix.  Propagation  of  the  known  n  x  1  matrix  state  estimate  X kk  is 
given  by: 

Xk+\\k  =  ^krl,k^k\k 

because  the  conditional  mean  of  Jk+l  k  is  zero.  Propagation  of  the  known  n  x  n  matrix  state  estimate  error 
covariance  matrix  Pk,k  is  given  by: 
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P  =0  P  0  +0 

1  k+\\k  ^k+\,kA  k\k^k+l,k  T  *Ik+\,k 

where  the  n  x  «  matrix  Qk+V[k  is  called  the  process  noise  covariance  matrix4  . 

Nonlinear  MU 

Calculate  the  n  x  1  matrix  filter  gain  Kk+i : 

K-k+ 1  =  Pk+\\k^k+\  \_H k+\Pk+\\kH k+\  +  ^<t+l] 

The  filter  measurement  update  state  estimate  n  x  1  matrix  Xk+l,k+x ,  due  to  the  observation  yk+l ,  is  calculated  with: 

TVu  —  T(^r+i|i) 


I— i 


^k+l\k+l  Xk+\\k  +  Kk+l 


where  Rk+l  is  the  scalar  variance  on  the  observation  residual  yk , ,  -y(xk+ and  y[xk+x ^  )  is  a  nonlinear 
function  of  Xk+l,k  .  Define  the  error  AXi+1,i+1  inXk+^k+l : 

^-k+\\k+\  =  Xk+\  —  ^i+l|i+l 
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Figure  6.  KF1  frequency  Errors. 


Define  the  n  x  n  state  estimate  error  covariance  matrix  Pk+V[k+ 1  with: 
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n»,|,«=£j(AXw|w)(AXM|w)r} 


Bucy  and  Joseph  [3]  (page  141)  recommend  that  Pk+uk+i  be  calculated  with: 


where: 

has  symmetry,  and: 


P  =P  -T 

1  k+\\k+\  1  k+\\k  1 


T  —  p  prT  d1  zj  p 

1  rk+\\k11  k+l^k+l11  k+lrk+\\k 


R,.„  =H,„,P„UH„,+R 


vi+ 1 


k+\*  k+\\k**  k+\  1  “A'+l 


Calculation  of  Pk+l ,t+1  by  the  last  three  equations  is  numerically  symmetric.  They  reduce  to  the  form  given  by 
Kalman: 

Pk+ 1|*+1  =  ~  ^k+\Hk+\  ]  Pk+\\k 

which  is  not  numerically  symmetric. 


KF1  Drift  Errors  (s/sA2) 


100000  300000  500000  700000 


time(s) 


Figure  7.  KF1  drift  errors. 


VI).  Nonlinear  TU  and  Nonlinear  MU 

Refer  to  Subsection  Nonlinear  MU  of  Section  Vc  for  the  nonlinear  MU. 
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Nonlinear  TU 


The  nonlinear  TU  always  spans  a  non-empty  time  interval  and  requires  the  use  of  a  numerical  state  estimate 
integrator  (p,. .  Given  an  initial  time  t0,  a  final  time  t  f  ,  and  a  force  model  u(x(t),t^,  then  (px  propagates  the 

state  estimate  X(/0)  from  t0  to  tf  using  forces  u[x{t),t )  to  get  x[t  That  is: 


X  (tf )  =  A  {*/ \X  (*o )  ’t0 ,  M  (x  ( r) ,  r )  ,t0  <  T  <  tf 

This  can  be  shortened  to  write: 

X  (tj-  )  =  (f>x  | tf ,  X  (Jq  ) ,  t0 1 

where  the  use  of  force  function  u[x{t),t |  is  tacitly  implied.  Thus,  (px  is  a  column  matrix  with  n  elements: 


<t>.= 


6 

T  xn 


VE.  Kalman  Filter  Advantage 

Severe  computational  problems  are  incurred  in  any  attempt  to  estimate  unobservable  states  using  iterated  batch 
least-squares  methods  or  iterated  maximum  likelihood  methods  for  navigation,  because  state-sized  inversions  of 
singular  matrices  are  always  required.  Here,  the  Kalman  filter  is  distinguished  in  that  estimates  of  unobservable 
states  can  be  created  and  used  without  matrix  inversion  problems,  because  the  Kalman  filter  MU  is  free  of  state¬ 
sized  matrix  inversions. 

By  design,  one  typically  estimates  observable  states.  But  the  Kalman  filter  enables  one  to  create  unobservable 
states.  The  USAF  chose  to  create  unobservable  GPS  clock  parameter  states  for  construction  of  GPS  Time. 


VI.  KALMAN  FILTERS  KF1  AND  KF2 

1  have  simulated  GPS  pseudo-range  measurements  for  two  GPS  ground-station  clocks  SI  and  S2,  and  for  two 
GPS  NAVSTAR  clocks  N1  and  N2.  Here,  I  set  simulated  measurement  time  granularity  to  30  s  for  the  set  of  all 
visible  link  intervals.  Visible  and  non-visible  intervals  are  clearly  evident  in  Figure  8.  I  set  the  scalar  root- 

variance  y/R  for  both  measurement  simulations  and  Kalman  filter  KF1  to  Jr  =  1  cm.  Typically  Jr  =1  m  for 
GPS  pseudo-range,  but  when  carrier-phase  measurements  are  processed  simultaneously  with  pseudo-range,  the 

root-variance  is  reduced  by  two  orders  of  magnitude.  So  the  use  of  J R  =  1  cm  enables  me  to  quantify  lower 
performance  bounds  for  the  simultaneous  processing  of  both  measurement  types. 

VIA.  Create  GPS  Clock  Ensemble 

Typically,  one  processes  measurements  with  a  Kalman  filter  to  derive  sequential  estimates  of  a  multidimensional 
observable  state.  Instead,  here  I  imitate  the  GPS  operational  procedure  and  process  simulated  GPS  pseudo-range 
measurements  with  KF1  to  create  a  sequence  of  unobservable  multidimensional  clock  state  estimates.  Clock  state 
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components  are  unobservable  from  GPS  pseudo-range  measurements.  See  Figure  2  for  an  example  of  an 
ensemble  of  estimated  unobservable  clock  phase  deviation  state  components  created  by  KF1. 

Sherman's  Theorem 

GPS  Time,  the  unobservable  GPS  clock-ensemble  mean  phase,  is  created  by  the  use  of  Sherman's  Theorem 
[10,15]  in  the  USAF  Kalman  filter  measurement  update  algorithm  on  GPS  range  measurements.  Satisfaction  of 
Sherman's  Theorem  guarantees  that  the  mean-squared  state  estimate  error  on  each  obsen>able  state  estimate 
component  is  minimized.  But  the  mean-squared  state  estimate  error  on  each  unobsen’able  state  estimate 
component  is  not  reduced.  Thus,  the  unobservable  clock  phase  deviation  state  estimate  component  common  to 
every  GPS  clock  is  isolated  by  application  of  Sherman's  Theorem.  An  ensemble  of  unobservable  state  estimate 
components  is,  thus,  created  by  Sherman's  Theorem  -  see  Figure  2  for  an  example. 


VIb.  Initial  Condition  Errors 

A  significant  result  emerges  due  to  the  modeling  of  Kalman  filter  (KF1)  initial  condition  errors  in  phase, 
frequency,  and  frequency-drift.  Initial  estimated  clock  phase  deviations  are  significantly  displaced  by  the  KF1 
initial  condition  errors  in  phase.  As  time  evolves,  estimated  clock  phase  deviation  magnitudes  diverge 
continuously  and  increasingly  when  referred  to  true  (simulated)  phase  deviations,  and  this  is  due  to  filter  initial 
condition  errors  in  frequency  and  frequency  drift.  See  Figures  3  and  2  for  an  example. 

VIC.  PARTITION  OF  KF1  ESTIMATION  ERRORS 

Subtract  estimated  clock  deviations  from  simulated  (true)  clock  deviations  to  define  and  quantify  Kalman  filter 
(KF1)  estimation  errors.  Adopt  Brown’s  additive  partition  of  KF1  estimation  errors  into  two  components.  I  refer 
to  the  first  component  as  the  Unobservable  Error  Common  to  each  Clock  (UECC),  and  to  the  second  component 
as  the  Observable5  Error  Independent  for  each  Clock  (OEIC).  Initially,  the  variances  on  the  UECC  and  OEIC  are 
identical.  On  processing  the  first  GPS  pseudo-range  measurements  with  KF1,  the  variances  on  both  fall  quickly. 
But  with  continued  measurement  processing,  the  variances  on  the  UECC  increase  without  bound,  while  the 
variances  on  the  OEIC  approach  zero  asymptotically. 

For  simulated  GPS  pseudo-range  data,  I  create  an  optimal  sequential  estimate  of  the  UECC  by  application  of  a 
second  Kalman  filter  KF2  to  pseudo-measurements  defined  by  the  phase  components  of  KF1  estimation  errors. 

Since  there  is  no  physical  process  noise6  on  the  UECC,  an  estimate  of  the  UECC  can  also  be  achieved  using  a 
batch  least-squares  estimation  algorithm  on  the  phase  components  of  KF1  estimation  errors  -  demonstrated 
previously  by  Greenhall  [6,7]. 

VID.  UNOBSERVABLE  ERROR  COMMON  TO  EACH  CLOCK 

There  are  at  least  four  techniques  to  estimate  the  UECC  when  simulating  GPS  pseudo-range  data.  First,  one 
could  take  the  sample  mean  of  KF1  estimation  errors  across  the  clock  ensemble  at  each  time  and  form  a  sample 
variance  about  the  mean;  this  would  yield  a  sequential  sampling  procedure,  but  where  each  mean  and  variance  is 
sequentially  unconnected.  Second,  one  can  employ  Ken  Brown's  Implicit  Ensemble  Mean  (IEM)  and  covariance; 
this  is  a  batch  procedure  requiring  an  inversion  of  the  KF1  covariance  matrix  followed  by  a  second  matrix 
inversion  of  the  modified  covariance  matrix  inverse;  this  is  not  a  sequential  procedure.  Third,  one  can  adopt  the 
new  procedure  by  Charles  Greenhall  [7]  wherein  KF1  phase  estimation  errors  are  treated  as  pseudo¬ 
measurements,  and  are  processed  by  a  batch  least-squares  estimator  to  obtain  optimal  batch  estimates  and 
covariance  matrices  for  the  UECC.  Fourth,  one  can  treat  the  KF1  phase  estimation  errors  as  pseudo- 
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measurements,  invoke  a  second  Kalman  filter  (KF2),  and  process  these  phase  pseudo-measurements  with  KF2  to 
obtain  optimal  sequential  estimates  and  variances  for  the  UECC.  I  have  been  successful  with  this  approach. 

Figure  4  presents  the  SI,  S2,  Nl,  N2  ensemble  of  KF1  phase  estimation  error  functions,  overlaid  with  the  KF2 
sequential  UECC  estimated  function  in  phase.  The  KF2  UECC  estimated  function  is  appropriate  only  for  KF1 
clock  phase  errors,  not  for  frequency  or  frequency  drift.  Figures  6  and  7  present  KF1  error  functions  in  frequency 
and  frequency  drift. 

VIE,  OBSERVABLE  ERROR  INDEPENDENT  FOR  EACH  CLOCK 

At  each  applicable  time,  subtract  the  estimate  of  the  KF2  phase  UECC  from  the  KF1  phase  deviation  estimate,  for 
each  particular  GPS  clock,  to  estimate  the  OEIC  for  that  clock.  Figure  8  presents  a  graph  of  the  phase  OEIC  for 
ground  station  clock  SI.  Intervals  of  KF1  range  measurement  processing  are  clearly  distinguished  from 
propagation  intervals  with  no  measurements.  During  measurement  processing,  the  observable  component  of  KF1 
estimation  error  is  contained  within  an  envelope  of  a  few  parts  of  a  nanosecond. 


SI  Observable  Component  Phase  Errors  (s) 


Case  13 


Figure  8.  SI  observable  component  of  phase  errors. 


Calculation  of  the  sequential  covariance  for  the  OEIC  requires  a  matrix  value  for  the  cross-covariance  between 
the  KF1  phase  deviation  estimation  error  and  the  UECC  estimation  error  at  each  time.  I  have  not  yet  been  able  to 
calculate  this  cross-covariance. 
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VII.  OBSERVABILITY 

I  have  defined  observability  in  terms  of  a  Kalman  filter  formulation,  and  1  have  proved  simple  theorems  related 
thereto.  My  definition  of  observability  is  different  than  Kalman’s  definition  and,  unlike  Kalman’s  definition,  is 
directly  applicable  to  covariance  matrices  derived  from  a  Kalman  filter. 

VIlA.  DEFINITION 

If  the  state  estimate  error  variance  of  a  particular  state  estimate  component  is  reduced  by  processing  an 
observation,  then  that  state  estimate  component  is  obsen’able  to  that  observation.  Otherwise,  that  state  estimate 
component  is  not  obsen’able  ( unobservable )  to  that  observation. 

VIlB.  THEOREM  1 

If  every  component  of  the  row  matrix  Hk+1  of  measurement- state  partial  derivatives  is  zero  at  time  tk+]  ,  then 
every  component  of  the  state  estimate  Xk+l  is  unobsen’able  at  time  tk+l  . 

Proof 

Hk+ j  =  0  implies  that  Pk+«k+x  =  Pk+l ,k .  Thus  none  of  the  variances  of  Pk , ,  k  are  reduced  due  to  processing  the 
observation  J4+1 .  Then  by  definition,  Xk+l  is  unobservable  in  every  component. 

VIIC.  THEOREM  2 

Given  values  for  scalars  Hk+l,  Pk , ,  k  >  0,  Rk+l  >  0  at  time  tk+l ,  and  given  that  Hk+1  T  0  ,  then  the  scalar  state 
estimate  Xk+l  is  observable  at  time  tk+1 . 

Proof 

The  obvious  inequality  Pk+l*kHk+l  +  Rk+l  >  Pk+i\kHL\  >  0  implies  that: 


Multiply  through  by  - 1 : 


Add  1: 


Multiply  through  by  Pk+l ,k 


Then: 


P  H2 

1>  W^+!  >Q 

Pk+l\k^k+l  +  Pk+ 1 


p  H 

{< - k+i\k  k+ 1 - <Q 

Pk+\\k^k+\  +  Pk+ 1 


0< 


1- 


P  H 

1  k+Vik11  k+ 1 


Pk+\\k^k+\  +  R-k+l 


<1 


0< 


1- 


p  H 

1  k+\\k11k+\ 


Pk+l\kH k+l  +  &k+l 


P  <P 

1  k+\\k  ^  1  k+\\k 


382 


39th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


Therefore: 


k+\\k+\ 


1  — 


p  H 

1  k+\\k11k+\ 


^k+l\k^k+l  +  ^k+ 1 


k+\\k 


o  <  p  <  p 

u  ^  1  k+l\k+\  v  1  k+\\k 


Thus,  the  variance  Pk+l ,k  is  reduced  due  to  processing  the  observation  yk+l .  Then  the  scalar  state  Xk+1  is 
observable  by  definition. 

VIId.  Theoretical  Foundation 

These  theorems  are  referred  to  expressions  given  by  Kalman  for  filter  gain  Kk+i  and  covariance  Pk+l ,Jt+1 . 

Kalman’s  expressions  are  derived  from  the  rigorous  theorem  chain  provided  by  Sherman,  Anderson,  and  Brunn- 
Minkowski  -  the  theoretical  foundation  is  deep. 

VIlE.  DETERMINE  OBSERVABILITY  DIRECTLY 

Given  an  optimal  sequential  estimator,  given  a  particular  collection  of  applicable  observations  (real  or  simulated), 
and  given  realistic  state  estimate  error  covariance  matrices  Pk+1,k  and  Pk+l ,i+1  at  each  time  tk+l ,  apply  the 

definition  of  observability  directly7  to  distinguish  between  observable  and  unobservable  state  elements.  An 
optimal  sequential  estimator  is  designed  to  eliminate  significant  aliasing  between  estimated  state  elements  and, 
thus,  enable  this  distinction. 


VIII.  UNOBSERVABLE  GPS  CLOCK  STATES 

GPS  Time  is  created  by  the  operational  USAF  Kalman  filter  by  processing  GPS  pseudo-range  observations.  GPS 
Time  is  the  mean  phase  of  an  ensemble  of  many  GPS  clocks,  and  yet  the  clock  phase  of  every  operational  GPS 
clock  is  unobservable  from  GPS  pseudo-range  observations,  as  demonstrated  below.  GPS  NAVSTAR  orbit 
parameters  are  observable  from  GPS  pseudo-range  observations.  The  USAF  Kalman  filter  simultaneously 
estimates  orbit  parameters  and  clock  parameters  from  GPS  pseudo-range  observations,  so  the  state  estimate  is 
partitioned  in  this  manner  into  a  subset  of  unobservable  clock  parameters  and  a  subset  of  observable  orbit 
parameters.  This  partition  is  performed  by  application  of  Sherman’s  Theorem  in  the  MU. 

VIIIa.  GPS  Pseudo-Range  Representation 

Let  tjh  denote  time8  of  radio  wave  transmission  for  the  h"’  NAVSTAR  clock,  and  let  tf  denote  time  of  radio 
wave  receipt  for  the  i,h  ground  station  clock.  Let  Ax}1'  and  Ax}'  denote  Kalman  filter  estimation  errors  in  clock 
phase  for  t}h  and  t}‘  .  Define  time  of  transmission  difference  t}  and  time  of  receipt  difference  t} : 


Thus: 


, Dh  _  fNh 

V  y  t  y 

fDi  _  fGi 
lR  ~  lR 


-oxT 


,Nh  _  .Dh  . 

v  y  I  L/  Ay 
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tR‘  =tR‘  +  SxR 


Gi 


Define  the  one-way  GPS  pseudo-range  measurement  pNhGi : 

Pma,=4«-C\  tf  >t“ 

Then: 

Pmci  =  c  (  \}r  +  Sxr  ]  -  [t™  +  SxT  ] ) 
=  c  ( [ff  -  tf  ]  +  [Sxf  -  8x. f  ]) 
where  c  is  speed  of  light  in  vaccua.  Define: 

A t  =  tf  -tf 
8thi  =  Scf  -  Scf 

Then: 

PNhGi  =  +  &  ) 

where  At  is  deterministic  and  8t  is  random. 

VIIlB.  PARTITION  OF  KALMAN  FILTER  ESTIMATION  ERRORS 


Let  xc  denote  the  phase  component  of  Kalman  filter  estimation  error  that  is  common  to  every  GPS  ensemble 
clock,  when  it  exists.  Define  phase  differences  x((''R  and  Xq't  with: 


XOR 

Nh 

XQt 


-x 


c 


for  ground  station  i  and  NAVSTAR  h  .  Then  Kalman  filter  estimation  errors  8xR  ,  i  e  { 1, 2, . . .  j  ,  for  ground 
station  clocks  and  SxRh ,  h  e  [l,  2,. . .} ,  for  NAVSTAR  clocks  have  the  additive  partition9: 


<T  Gi  ,  Gi 

OXr  —  xc  +  xOR 

Sxnt"  =xc+x“ 


Vine.  The  Common  Random  Phase  Component  Is  Unobservable 

With  substitutions: 

8thi  =  Sxf  - 
=  [xc  +Xqr\~  [xc  +  Xqj  J 

_  Gi  _  Nh 
—  XQR  Xot 

Then: 
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PNhGi=c(At  +  \xOR  -*ori 

Thus,  the  random  phase  component  xc  that  is  common  to  the  Kalman  filter  estimation  error  for  every  ensemble 
clock  has  vanished  in  the  range  representation  pNhGi .  Variations  Axc  in  xc  cannot  cause  variations  A pNhGi  in 

Pma  • 

An  -  df>mGi  Ax 
A PNhGi  t-SXc 

dxc 

because  the  partial  derivative  H  =  dpmGi  /  Gxc  is  zero: 

^PmGi  _  Q 

5xc 


I  have  thus  shown  that  xc  is  unobservable  from  pNhGi .  But  the  architect  who  designs  the  complete  estimator 
must  design  an  optimal  NAVSTAR  orbit  estimator  to  prevent  aliasing  from  NAVSTAR  orbit  estimation  errors 
into  xc  .  It  helps  to  know  that  there  is  no  coupling  between  the  orbit  and  xc  in  the  complete  state  transition 

function.  I  have  provided  a  new  method  herein  to  identify  this  aliasing,  and  I  have  provided  suggestions  on  where 
to  look  for  inadequate  modeling  that  would  be  the  source  of  this  aliasing.  See  Section  X. 

VIIlD.  INDEPENDENT  RANDOM  PHASE  COMPONENTS  ARE  OBSERVABLE 


The  independent  phase  deviations  xGR  and  xG‘r  are  observable  to  GPS  pseudo-range  observations  because  their 
partial  derivatives  are  non-zero: 

^PmlL  =  +c 
dx?1 


v  OR 


op 


NhGi 


dx™ 


=  -c 


Estimation  of  xGR  and  xG'T  by  the  Kalman  filter  will  reduce  their  error  variances. 


IX.  ALLAN  VARIANCE  AND  PPN  RELATIONS 

IXa.  Allan  Coefficients  vs.  Diffusion  Coefficients  For  2-State  Clocks 

Denote  r  as  clock  averaging  time,  o\  ( r  j  as  Allan  variance,  a0  as  Allan’s  FMWN  coefficient,  a  2  as  Allan’s 
FMRW  coefficient,  <r]  as  the  FMWN  diffusion  coefficient,  and  c1  as  the  FMRW  diffusion  coefficient.  Then: 


where: 


2  (  \  -1  2-1 
Gy  [T )  =  a0T  +  a_2T  =  G{T 


1  2 

+— a\r 
3 
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c 2  —  -J3a_2 

Notice  that  a\  (r)  depends  only  on  the  time  difference  z  =  tk+l  —  tk ,  not  on  the  time  tk  for  2-state  clocks. 


Proportionate  Process  Noise  (PPN) 

Let  a  denote  a  variable  a  e  [l,  2, 3, . . . ,  TV  j  to  identify  each  GPS  clock  in  an  ensemble  of  N  clocks.  For  each 
clock  a  ,  define  the  ratio  Sa  between  diffusion  coefficients  crla  and  o2a : 


cr, 


la 


Then  PPN  is  defined  when,  for  each  GPS  clock  a  and  each  associated  ratio  Sa ,  we  have: 


sl  =  s2=s3  =  -  =  sN 


IXb.  Allan  Coefficients  vs.  Diffusion  Coefficients  For  3-State  Clocks 

The  reader  is  referred  to  Zucca-Tavella  [18]  Equation  43  for  the  time-dependent  representation  of  the  Allan 
variance: 


i/  \  2-1  It  1  23 

;(r)  =  07 r  ‘+-cr;r  + — at 
w  1  3  20  3 


+  cr; 


1  3  1  2  "i 

—  T  +—Tt, 

3  2  k) 


+  2r2(C3  +^3[r  +  ^])2 


where  c3  and  ju3  denote  initial  condition  and  deterministic  mean  for  frequency  drift. 

Case  13 

For  Case  13,  the  diffusion  coefficient  values  crla  and  ola  are  same  for  each  simulated  clock  Sa  =  <J2a  /  crla , 
a  e  |l,2, 3, 4}  : 

S  =1.667  xlO-4  s  1 


But  PPN  has  not  been  defined  for  3-state  clocks. 


X.  IDENTIFY  NON-CLOCK  MODELING  ERRORS 

My  interest  in  the  GPS  NAVSTAR  (SV)  orbit  determination  problem,  combined  with  that  of  the  clock  parameter 
estimation  problem,  has  enabled  the  identification  of  a  useful  diagnostic  tool:  Given  realistic  values  for  diffusion 
coefficients  for  each  of  the  real  GPS  clocks,  then  quantitative  upper  bounds  can  be  calculated  on  OEIC 
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magnitudes.  These  calculations  require  the  use  of  a  rigorous  simulator.  Existence  of  significant  cross¬ 
correlations  between  GPS  clock  phase  errors  and  other  non-clock  GPS  estimation  modeling  errors  enables 
significant  aliasing  into  GPS  clock  phase  estimates  during  operation  of  KF1  on  real  data.  But  given  rigorous 
quantitative  upper  bounds  on  OEIC  magnitudes,  then  significant  violation  of  these  bounds  when  processing  real 
GPS  pseudo-range  and  carrier-phase  data  identifies  non-clock  modeling  errors  related  to  the  GPS  estimation 
model.  Modeling  error  candidates  here  include  NAVSTAR  orbit  force  modeling  errors,  ground  antenna  modeling 
errors  (multipath),  and  tropospheric  modeling  errors.  NAVSTAR  orbit  force  modeling  errors  include  those  of 
solar  photon  pressure,  albedo,  thermal  dump,  and  propellant  outgassing.  The  accuracy  of  this  diagnostic  tool 
depends  on  the  use  of  realistic  clock  diffusion  coefficient  values  and  a  rigorous  clock  model  simulation  capability. 


XI.  UNOBSERVABLE  MEAN  PHASE 

In  an  earlier  version  of  my  paper,  I  reported  on  KF1  validation  results  where  clock  SI  was  specified  as  a 
TAI/UTC  clock,  external  to  the  GPS  clock  ensemble  consisting  of  S2,  Nl,  and  N2.  This  brought  observability  to 
S2,  Nl,  and  N2  clock  states  from  GPS  pseudo-range  measurements,  drove  clocks  S2,  Nl,  and  N2  immediately  to 
the  TAI/UTC  timescale,  and  enabled  a  clean  validation  of  my  filter  implementation.  Also,  it  raised  the  question: 
Why  not  do  the  same  thing  for  the  real  GPS  clock  ensemble?  Discussions  with  Ed  Powers  (USNO)  and  Bill 
Feess  (Aerospace  Corporation)  reveal  that  this  approach  was  tried  and  discarded  after  the  difficulty  in  recovery 
from  an  uplink  hardware  failure  was  blamed  on  the  use  of  a  single  TAI/UTC  Master  Clock.  This  issue  was 
resolved  with  Kenneth  Brown’s  introduction  of  the  implicit  ensemble  mean.  The  mean  phase  (GPS  Time)  of  the 
GPS  clock  ensemble  will  remain  unobservable  to  GPS  pseudo-range  measurements  in  the  USAF  Kalman  filter  for 
the  foreseeable  future. 


NOTES 

1.  James  R  Wright  is  the  architect  of  ODTK  (Orbit  Determination  Tool  Kit),  a  commercial  software  product 
offered  by  Analytical  Graphics,  Inc.  (AGI). 

2.  Session  IX,  Paper  34. 

3.  According  to  Bill  Feess,  an  improvement  in  control  can  be  achieved  by  replacing  the  existing  “bang-bang 
controller”  with  a  “proportional  controller.” 

4.  See  Zucca  and  Tavella  [18]  for  concrete  clock  examples  of  Jk+i  k  and  Qk+t  k  . 

5.  Observability  is  meaningful  here  only  when  processing  simulated  GPS  pseudo-range  data. 

6.  I  apply  sufficient  process  noise  covariance  for  KF2  to  mask  the  effects  of  double-precision  computer 
word  truncation.  Without  this,  KF2  does  diverge. 

7.  Note  that  this  is  impossible  using  Kalman’s  definition  of  observability. 

8.  Refer  all  times  to  a  coordinate  time,  e.g.,  to  GPS  Time.  Appropriate  transformations  between  proper  time 
and  coordinate  time  must  be  performed  in  the  operational  algorithms,  but  state  estimate  observability  is 
independent  of  relativity,  so  observability  can  be  defined  and  discussed  independent  of  relativity. 

9.  This  partition  was  introduced  by  Kenneth  Brown  [2]. 
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